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Large-eddy simulation of flow ~ 

in a plane, asymmetric diffuser 

By Hans-Jakob Kaltenbach 

Recent improvements in subgrid- scale modeling as well as increases in computer 
power make it feasible to investigate flows using large-eddy simulation (LES) which 
have been traditionally studied with techniques based on Reynolds averaging. How- 
ever, LES has not yet been applied to many flows of immediate technical interest. 
This report describes preliminary results from LES of a plane diffuser flow. The 
long term goal of this work is to investigate flow separation as well as separation 
control in ducts and ramp-like geometries. 

1« Motivation and objectives 

Flow separation is a fundamental phenomenon which occurs in many engineering 
flow configurations such as airfoils, diffusers, and cascades. Prevention of separation 
usually improves the performance of such devices by increasing pressure recovery, 
enhancing lift, and decreasing total drag. Methods which delay or prevent sepa- 
ration include passive devices like vortex generators and active systems which use 
auxiliary power to modify the flow. Recent experiments (Katz et al 1989 , Seifert 
el al. 1992, 1993, Obi et al. 1993) have shown that flow control strategies which 
rely on adding small amounts of periodic disturbances (e.g. suction and blowing at 
specific locations) can delay separation efficiently without consuming large amounts 
of auxiliary power. 

Investigation of unsteady control concepts using numerical simulation requires 
a method which computes the spatial as well as the temporal evolution of the 
turbulent fluctuations. Methods based on the Reynolds averaged Navier Stokes 
equations have difficulty dealing with the unsteadiness which is essential for the 
control concept. LES seems to be an ideal tool for studying unsteady control be- 
cause experimental work has shown that the delay of separation is mostly due to the 
creation of large-scale vortical structures which improve the entrainment of high- 
speed fluid into separated zones (Katz et al. 1989). LES is very likely to adequately 
describe this physical mechanism at a reasonable cost. Conversely, direct numer- 
ical simulation (DNS) would be lather expensive at the Reynolds numbers under 
consideration. 

The use of generalized coordinates largely enhances the range of possible flow 
configurations accessible for LES. Finite differences are more convenient to use 
than highly accurate spectral methods for approximation of derivatives. Little is 
known about the mutual dependency of use of the SGS-model in combination with 
a numerical scheme with considerable dispersive properties like second order cen- 
tral differences. One way of exploring this dependency consists of studying the 
sensitivity of the solution with respect to grid resolution. 
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FIGURE 1. Computational domain for the plane diffuser in units of 6. Only each 
4th of the wall-normal and each third of the streamwise grid lines are plotted. 

The objective of this study is to perform large-eddy simulation of turbulent flow 
in a one-sided diffuser. The results will be validated by comparison with the experi- 
ment of Obi et al (1993). Once the minimum resolution requirements for accurately 
simulating this flow have been determined, periodic disturbances with various fre- 
quencies and amplitudes will be added to the flow as was done in the experiment. 
The final goal of this study is to understand the physics of flow control through 
periodic forcing and to determine the parameters (location, frequency, amplitude) 
that are most efficient in delaying flow separation. 

2. Accomplishments 

The DNS code of Choi & Moin (1993) has been modified to allow for simulation 
of inflow/outflow configurations with either no-slip or no-stress conditions at the 
upper boundary. A detailed description of the coordinate transformation and the 
numerical scheme is given in Choi et al. (1992). The code solves the unsteady, 
incompressible Navier-Stokes equations in generalized coordinates in two dimensions 
and a Cartesian equidistant grid in the third (spanwise) direction. The fully-implicit 
time integration scheme (Crank- Nicholson) uses Newton linearization along with 
approximate factorization. The time advancement allows rather large time-steps 
and CFL- numbers. The Poisson solver makes use of the spanwise periodicity of the 
domain which allows a Fourier transformation in this direction. The remaining 2D 
problems are solved by using LU decomposition for the first spanwise wavenumber 
and an iterative solver for higher wavenumbers. The major part of CPU-time needed 
for advancing the fractional step scheme by one timestep is spent on performing 3-4 
Newton iterations per timestep. Depending on the spanwise resolution, the Poisson 
solver consumes between 10 and 30% of the computational effort. 

Unsteady data are specified at the inflow, and a convective boundary condition is 
applied at the outflow. The upper and lower boundaries are no-slip walls. The data 
to be specified at the inflow plane of the domain are created in an independent LES 
of a fully developed channel flow. At each timestep, a cross-section of the velocity 
field is stored on disk using data reduction to 4 Byte words in order to minimize 
data storage and input/output time. These channel flow data are subsequently fed 
into the code which simulates flow through a non-periodic duct configuration. 

A simple version of the dynamic SGS model (Germano et al. 1991) has been 
implemented and tested for a periodic channel flow. It makes use of least-square 
contraction (Lilly 1992) in combination with spanwise averaging. The total viscosity 
is constrained to be positive through a clipping operation. The test filter is applied 
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in the spanwise direction and along horizontal lines in the transformed space. No 
filtering takes place in the wall-normal direction. The SGS-model increases the 
CPU-time by less than 10% and has no significant effect on memory requirements. 
A formulation of the SGS-model in generalized coordinates is given in the appendix. 

The code has been used for simulation of turbulent flow through a plane, one- 
sided diffuser. The dimensions of the computational domain are shown in figure 
1. The diffuser geometry and the Reynolds number Ret, = Ub8/u = 9000 match 
the experimental configuration of Obi et ai (1993). Here, Ub denotes the bulk 
velocity of the incoming fully developed turbulent channel flow of height 28. The 
flow from the inlet channel of length 6(5 enters an asymmetric diffuser with an 
expansion ratio of a = 4.7. The diffuser and the outlet channel extend over 42 8 
and 52 6, respectively. The upper wall remains parallel to the inlet channel whereas 
the lower wall is deflected by approximately 10°. Both corners formed by the inlet 
and outlet channels with the lower wall are slightly rounded. With a width of 12 6, 
the aspect ratios of inlet and outlet channel are 1 : 6 and 1 : 1.28, respectively. 
The experiment had much higher aspect ratios of 1 : 35 and 1 : 7.45 in order to 
guarantee a two-dimensional core flow free from sidewall effects. 

The grid spacing is based on an estimate for the skin friction coefficient as a 
function of bulk Reynolds number as given by Dean; cj t b = 0.061 i?e^ 0,25 . For 
Reb = 9000, this relation gives a value c/ b = 0.0063, which corresponds to Re T = 
u r 6/ v = 500. Therefore, a wall- unit in the inlet section of the diffuser is 0.002(5. 
Assuming that the flow in the outlet channel will finally evolve into a fully developed 
channel flow, the wall-unit in the outflow section increases by the expansion ratio, 
a, to 0.0094(5. Linear stretching of the mesh sizes in the streamwise and wall-normal 
directions accounts for this increase. The grid with Ax = 0.375 8 (Ax + = 190) in 
the inflow and Ax = 1.25 8 (Ax + = 133) in the outflow section has 124 streamwise 
points. In the wall-normal direction, 64 points are distributed by using hyberbolic 
tangent stretching. It should be noted, however, that estimates for grid spacing in 
the outflow section based on wall- units are probably not very relevant because the 
flow is partially separated. The proper spacing must be verified from sensitivity 
studies. 

For the present flow it would be highly desirable to use “zonal grids” for the 
spanwise direction. This would allow accurate simulation of the inlet with a fine 
spanwise grid and with increasingly coarser grids towards the diffuser outlet. How- 
ever, the spanwise spacing A z has to remain constant throughout the domain in 
the present code. The spanwise spacing must be chosen as a compromise between 
the different requirements in the inflow and the outflow section. Distribution of 64 
points over a spanwise extent of 12 8 results in A z+ = 94 in the inflow and Ar + = 20 
in the outflow section. These estimates are based on Re r = 500, which might not 
be reached everywhere in the actual simulation. The grid is fine enough to resolve 
near wall streaks with a characteristic spacing of approximately 100 wall-units in 
the diffuser outlet but certainly too coarse to realistically simulate the flow in the 
diffuser inlet section. For the purpose of studying the effect of unsteady blowing 
and suction on the pressure recovery in the rear part of the diffuser, the features 
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FIGURE 2. (a) Pressure coefficient c p along the upper wall and (b) skin friction 

c / based on centerline velocity at the inflow plane along lower ( ) and upper 

( ) walls as function of streamwise distance from the inflow plane. 

of the flow in the inlet channel might be of minor importance. However, this issue 
must be tested by means of grid refinement studies. 

An estimate of the quality of the simulated inlet flow can be obtained by com- 
parison with Dean’s empirical relations for skin friction and ratio of the centerline 
to the bulk velocity. Dean’s correlation U c /Ub = 1.27iJe^ 0,0116 gives 1.14 for the 
Reynolds number considered. Three channel flow LES have been carried out for 
Ret = 9000 with spanwise spacings A = 12, 29 and 73 based on the actual wall 
shear which changes significantly with the grid resolution. Only the finest case 
reaches a Cf^ = 0.0064 (iZe r = 510), which is close to the value derived from 
Dean’s relation. The other cases reach c/ t b = 0.0054 (Re T = 465) and 0.0038 
( Re r = 393). This means that the coarsest simulation which is used as database 
for unsteady inflow underpredicts the skin friction by 40%. However, the values of 
Uc/Ub = 1.12, 1.09 and 1.10 do not differ much between the three cases, indicating 
that the interior of the flow is well captured by all three simulations. The case with 
the finest spacing is actually not too far from a DNS because the SGS-eddy viscos- 
ity contributes less than 20% to the total viscosity. The plane averaged SGS-eddy 
viscosity is larger than the molecular viscosity in 60% of the domain in the coarsest 
case. 

The inertial time scale r = Q.5h(x)/Ub(x) increases with the square of the expan- 
sion ratio from the inlet to the outlet section, i.e. r out = a 2 r in . Here, h(x) is the 
diffuser height at location x. The same holds for the viscous time scale r v i 3C = v/u 2 r 
which mainly determines the choice of the maximum timestep. Again - as for the 
spanwise spacing - the choice of the timestep is a compromise between the very 
different requirements of inlet and outlet section. Running the code with 10 time 
steps per Tj n corresponds to a timestep which is three times bigger than the viscous 
timescale of the inflow but only one seventh of the viscous timescale of the outflow. 
In order to obtain converged statistics, data have to be sampled over a period of 
100 inertial time scales. Because the inertial timescale of the flow increases by a 
factor of 22 from inlet to outlet rather long simulation times - of the order of 10000 
timesteps - are required. 
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FIGURE 3. Region with reversed flow. Shown are isolines of the contravariant 
velocity component which is aligned with the horizontal grid lines. The grid lines 
shown are a subset of the actual grid. 


Preliminary results from a coarse grid simulation are shown in figures 2 to 5. 
Data have been sampled over 300 inlet inertial timescales after running the code 
for 1000 inlet timescales. Statistics are not fully converged, but the main features 
of the flow seem to be established. 

Figure 2 depicts pressure recovery along the upper wall and skin friction along 
both walls. Negative values of Cf behind the inlet corner mark the streamwise 
extent of a first separation bubble. This separated zone, which is shown in more 
detail in figure 3, was not reported by Obi et ai (1993). The flow separation close 
to the inlet might be an artifact of the coarse spanwise resolution in the diffuser 
inlet. This can be understood in the following way: the low value of skin friction 
in the inlet channel indicates that the velocity profile is less full than expected for 
the given Reynolds number based on centerline velocity. A less full mean velocity 
profile arises because the turbulent shear stress is underpredicted. Less high speed 
fluid from the channel core is transported towards the wall which makes separation 
more likely to occur. Grid refinement studies will be carried out to clarify this point. 
The skin friction approaches zero in the vicinity of the diffuser outlet. In this region 
the flow seems to be close to separation. However, we do not find a zone of reversed 
flow as for the first bubble. In accordance with this finding, mean profiles in the 
experiment do not show significant backflow at the diffuser outlet. 

Profiles of mean velocity and turbulence intensities change drastically inside the 
diffuser. As described by Simpson (1985), relative high turbulence levels might be 
found in regions with backflow. Figure 4 shows a strong increase of the stream- 
wise velocity fluctuation in the vicinity of the first separation bubble (x/S « 9). 
Throughout the rear part of the diffuser, the turbulence intensities remain high 
compared to the value of the mean velocity (figure 5). The scatter in these pro- 
files is higher than in figure 4 because much longer sampling times are required in 
the outflow section to obtain converged statistics. There is little evidence from the 
mean flow profiles that the flow is separated at the diffuser exit. 
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FIGURE 4. Profiles of mean velocity U /U c ( ) and contravariant fluctuations 

v! fUc ( ), v * /Uc ( ) and w 1 jU c ( ) close to the diffuser inlet. Rms~ 

values are enhanced by a factor of five compared to the mean velocity. 
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FIGURE 5. Same as figure 4 for the rear part of the diffuser. Note the scale change 
between this figure and figure 4. 
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3. Future plans 

Grid refinement studies will be carried out in order to see whether underprediction 
of skin friction in the inlet channel has a significant effect on the flow inside the 
diffuser. A detailed comparison with available experimental data as well as results 
from simulations based on statistical turbulence models will be carried out. Once 
the undisturbed flow is predicted well, oscillatory control of the diffuser will be 
studied in order to see whether the same gains in pressure recovery are reached in 
the simulation as in the experiment. 

In addition to flow through a diffuser, flow separation at a corner formed by a 
wall with a hinged rearward facing ramp will be investigated. This configuration 
resembles the generic flap which is the subject of an ongoing experimental investi- 
gation (Seifert et al. 1993). The computational domain will be designed such that 
it matches the experimental configuration of Katz et al (1989). The long term goal 
of these studies, which are carried out in collaboration with Dr. I. Wygnanski and 
Dr. A. Seifert, Tel Aviv University, is the improvement of high-lift devices such as 
highly deflected flaps. 

Appendix: The dynamic SGS-model in generalized coordinates 

The dynamic model for residual stresses (Germano et a/., 1991) requires compu- 
tation of the Leonard term and the deformation rate tensor S ,Jf . One obtains 
the scalar model coefficient through contractions of these tensors (and their filtered 
counterparts). Several possibilities exist to formulate the model in generalized co- 
ordinates. If the model formulation is based on a Cartesian system, the Cartesian 
tensors L J and S lJ have to be computed from the contravariant velocity compo- 
U { jNj P roce< ^ ure ^ to moderately complex expressions for the components 

On the other hand, the model can be formulated in generalized coordinates. In 
this case, contravariant as well as covariant components of strain-rate and Leonard 
term A* J , A l; have to be evaluated. It turns out that this formulation with 

certain assumptions leads to rather simple expressions for the Leonard term. This 
will be shown in the following sections where the definitions 
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Cartesian coordinate 
generalized coordinate 
Cartesian velocity component 
contravariant velocity component 

velocity component weighted with Jacobian (mesh volume) 
aerivative of Cartesian with respect to generalized coordinates 
derivative of generalized with respect to Cartesian coordinates 
contravariant metric coefficient weighted with Jacobian 


are used. 
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Starting from the Cartesian formulation for the Leonard term L' } = v'v) - v'v>, 
where x stands for filtering, the contravariant tensor can be written as 

A «j _ (cj n ) -l (c“) -1 £ mn = (0 - H c ") -1 ( v "‘« n “ v™v") 

Under the assumption that cj changes only weakly over the distance corresponding 
to the test filter width, it might be extracted from beneath the filter operator. This 

leads to 

A ,J w (c™)- 1 (c") _1 c^c"(ui'u« - uPui) 

= S'pS^xiP u* - W) = (u’u J - u’u-j) . 

Thus, the Leonard term is computed in generalized coordinates from the contravari- 
ant velocity components in the same way as L' } from the Cartesian velocity com- 
ponents. 

For computation of the strain-rate, the use of non-conservative expressions for 
derivatives leads to compact expressions. For the purpose of computing the SGS- 
eddy viscosity, it is probably of little importance whether or not the derivatives 
are formulated strictly conservative. The code conserves momentum as long as the 
momentum fluxes over cell surfaces sire formulated with a “telescoping scheme. 
One obtains the contravariant components of the strain rate tensor from 


s” = (O' 


J ' 2 dy" dy m ' 


=(cr)-(c;)-i|(c;)-^ + W)-|^) 


“ T7h c < > a 8*v 8xi h 


2J 1 


dxP 


dxi 


where . 

v m = -jC™q l and v n = — 

J ** 

has to be inserted because the code is actually formulated in terms of the weighted 
velocities q'. 
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